Ειδικά Μαθήματα Βοτανικής
Ειδικά Μαθήματα Βοτανικής
- 1 Προετοιμασία
- 2 Δεδομένα προς ανάλυση
- 3 Προκαταρκτικές αναλύσεις
- 4 Κυρίως ανάλυση
- 4.1 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο enter
- 4.2 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο της βηματικής πολλαπλής γραμμικής παλινδρόμησης
- 4.3 Ποιό μοντέλο είναι καλύτερο;
- 4.4 Δημιουργία μιας λειτουργίας για τον υπολογισμό της σχετικής σημαντικότητας των μεταβλητών
- 4.5 Υπολογισμός της σχετικής σημαντικότητας των μεταβλητών
- 4.6 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο των ολικών υποσυνόλων
- 4.7 Απεικόνιση των αποτελεσμάτων του προηγούμενου βηματος
- 4.8 Έλεγχος των αποτελεσμάτων της παλινδρόμησης
- 4.9 Έλεγχος πιθανοτήτων
- 5 Εργασία για το σπίτι
Όλα τα απαιτούμενα δεδομένα και αρχεία θα τα βρείτε στον ιστότοπο του μαθήματος στο e-class
1 Προετοιμασία
1.1 Εγκατάσταση βιβλιοθήκης
Εγκατάσταση του πρώτου και κύριου πακέτου (και ενός βοηθητικού)
(Αυτό γίνεται μια φορά σε κάθε Η/Υ - εάν έχει εγκατασταθεί δεν χρειάζεται να δοθεί ξανά η εντολή αυτή)
## CRAN repository:
install.packages("easypackages", repos = "http://cran.us.r-project.org")
## Github:
devtools::install_github("drsimonj/twidlr")1.2 Φόρτωση βιβλιοθήκης
Φόρτωση του εν λόγω πακέτου
library(easypackages)1.3 Εγκατάσταση βιβλιοθηκών
Με την εντολή packages('insert_package_name_here') γίνεται η εγκατάσταση των πακέτων που μας ενδιαφέρουν1
Δοκιμάστε το
1.4 Φόρτωση βιβλιοθηκών
Εντολή για φόρτωση των πακέτων
libraries("car", "psych", "tidyverse", "broom", "twidlr", "leaps", "relaimpo",
"gvlma", "MASS")1.5 Working Directory
Με την ακόλουθη εντολή, βλέπουμε σε ποιόν φάκελο δουλεύουμε:
getwd()1.6 Αλλαγή του working directory
Ενώ με την εντολή αυτή, αλλάζουμε τον φάκελο στον οποίο δουλεύουμε και αποθηκεύουμε δεδομένα
[π.χ., setwd("E:/") εάν θέλουμε να δουλεύουμε στον σκληρό δίσκο Ε]
setwd("E:/Academic/Lessons/<c5><e9><e4><e9><ea><dc> <cc><e1><e8><de><ec><e1><f4><e1> <c2><ef><f4><e1><ed><e9><ea><de><f2>/Labs/Labs")2 Δεδομένα προς ανάλυση
2.1 Εισαγωγή αρχείου
Τώρα ας εισάγουμε το αρχείο το οποίο περιέχει τα δεδομένα τα οποία μας ενδιαφέρουν να αναλύσουμε
Aegean <- readxl::read_excel("G:/Academic/Lessons/EMB/Labs/Labs/Aegean Islands.xlsx")Κάθε αρχείο πρέπει να έχει ένα όνομα στην R
Ποιό είναι το όνομα του αρχείου μας;
2.2 Δομή των δεδομένων μας
Ας δούμε τα δεδομένα τα οποία εισάγαμε στην R
Aegeanκαι τι είδους δεδομένα είναι αυτά (ποιοτικά; ποσοτικά;)
str(Aegean)## Classes 'tbl_df', 'tbl' and 'data.frame': 59 obs. of 14 variables:
## $ Island : chr "Alonissos" "Amorgos" "Anafi" "Andros" ...
## $ Area : num 65 121 38 380 19.7 ...
## $ Elevation: num 475 823 584 1003 378 ...
## $ MAT : num 111 133 139 130 139 ...
## $ MAP : num 489 498 498 393 362 ...
## $ Temp : num 161 175 181 170 181 ...
## $ Rain : num 495 528 523 463 402 ...
## $ Dis : num 4.76 10.33 21.14 9.29 1.7 ...
## $ Dm : num 41 117.8 207.1 55.6 112.2 ...
## $ Geology : num 5 5 8 6 3 1 5 4 1 10 ...
## $ Total : num 515 728 504 915 311 ...
## $ Regional : num 21 75 45 58 28 8 24 52 34 96 ...
## $ Endemics : num 18 59 35 47 25 6 15 36 17 30 ...
## $ SIE : num 0 1 0 3 1 0 0 1 0 1 ...
Μπορούσαμε να πάρουμε την πληροφορία αυτή με άλλον τρόπο;
2.3 Λιγότερα δεκαδικά
Με την ακόλουθη εντολή, θα εμφανίζονται πλέον λιγότερα δεκαδικά:
options(digits = 2)3 Προκαταρκτικές αναλύσεις
3.1 Έλεγχος κανονικότητας
3.1.1 Α
Δημιουργούμε ένα αρχείο το οποίο ελέγχει εάν όλες οι μεταβλητές μας (εκτός της πρώτης στήλης, η οποία είναι ποσοτική μεταβλητή - το όνομα των νησιών), έχουν κανονική κατανομή:
normality <- lapply(Aegean[, 2:14], shapiro.test)3.1.2 Β
Ας δούμε τα αποτελέσματα του προηγούμενου βήματος
normality## $Area
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.3, p-value = 3e-15
##
##
## $Elevation
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.8, p-value = 2e-06
##
##
## $MAT
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.8, p-value = 1e-07
##
##
## $MAP
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.9, p-value = 2e-04
##
##
## $Temp
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.8, p-value = 3e-07
##
##
## $Rain
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.9, p-value = 5e-05
##
##
## $Dis
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.8, p-value = 5e-07
##
##
## $Dm
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.9, p-value = 0.01
##
##
## $Geology
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.9, p-value = 1e-04
##
##
## $Total
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.9, p-value = 3e-04
##
##
## $Regional
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.6, p-value = 2e-11
##
##
## $Endemics
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.5, p-value = 3e-13
##
##
## $SIE
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.2, p-value = 3e-16
Οι μεταβλητές με πιθανότητα μικρότερη από 0.05 ΔΕΝ εμφανίζουν κανονική κατανομή.
3.2 Έλεγχος συσχέτισης
Ας ελέγξουμε τη συσχέτιση μεταξύ των μεταβλητών μας
Εάν κάποιες έχουν συσχέτιση \(\geq 0.8-0.9\), πρέπει να τις αφαιρέσουμε από τη μεταγενέστερη ανάλυση (Γιατί;)
corr.test(Aegean[, 2:14], method = "spearman")## Call:corr.test(x = Aegean[, 2:14], method = "spearman")
## Correlation matrix
## Area Elevation MAT MAP Temp Rain Dis Dm Geology
## Area 1.00 0.75 -0.39 0.10 -0.42 0.07 0.15 -0.41 0.79
## Elevation 0.75 1.00 -0.33 0.07 -0.38 0.04 0.04 -0.32 0.67
## MAT -0.39 -0.33 1.00 0.00 0.96 0.05 -0.23 0.46 -0.32
## MAP 0.10 0.07 0.00 1.00 0.05 0.99 0.40 -0.26 0.03
## Temp -0.42 -0.38 0.96 0.05 1.00 0.09 -0.19 0.33 -0.36
## Rain 0.07 0.04 0.05 0.99 0.09 1.00 0.38 -0.22 0.01
## Dis 0.15 0.04 -0.23 0.40 -0.19 0.38 1.00 -0.27 0.22
## Dm -0.41 -0.32 0.46 -0.26 0.33 -0.22 -0.27 1.00 -0.23
## Geology 0.79 0.67 -0.32 0.03 -0.36 0.01 0.22 -0.23 1.00
## Total 0.95 0.77 -0.42 0.09 -0.43 0.05 0.16 -0.47 0.79
## Regional 0.63 0.65 -0.19 0.11 -0.24 0.10 0.10 -0.17 0.58
## Endemics 0.43 0.45 -0.03 -0.15 -0.12 -0.15 -0.07 0.18 0.45
## SIE 0.62 0.53 -0.31 0.11 -0.32 0.09 0.35 -0.32 0.50
## Total Regional Endemics SIE
## Area 0.95 0.63 0.43 0.62
## Elevation 0.77 0.65 0.45 0.53
## MAT -0.42 -0.19 -0.03 -0.31
## MAP 0.09 0.11 -0.15 0.11
## Temp -0.43 -0.24 -0.12 -0.32
## Rain 0.05 0.10 -0.15 0.09
## Dis 0.16 0.10 -0.07 0.35
## Dm -0.47 -0.17 0.18 -0.32
## Geology 0.79 0.58 0.45 0.50
## Total 1.00 0.67 0.48 0.67
## Regional 0.67 1.00 0.77 0.53
## Endemics 0.48 0.77 1.00 0.49
## SIE 0.67 0.53 0.49 1.00
## Sample Size
## [1] 59
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## Area Elevation MAT MAP Temp Rain Dis Dm Geology Total
## Area 0.00 0.00 0.11 1.00 0.06 1.00 1.00 0.06 0.00 0.00
## Elevation 0.00 0.00 0.47 1.00 0.14 1.00 1.00 0.53 0.00 0.00
## MAT 0.00 0.01 0.00 1.00 0.00 1.00 1.00 0.02 0.54 0.06
## MAP 0.46 0.61 0.98 0.00 1.00 0.00 0.08 1.00 1.00 1.00
## Temp 0.00 0.00 0.00 0.68 0.00 1.00 1.00 0.47 0.25 0.04
## Rain 0.58 0.77 0.73 0.00 0.50 0.00 0.14 1.00 1.00 1.00
## Dis 0.27 0.75 0.08 0.00 0.15 0.00 0.00 1.00 1.00 1.00
## Dm 0.00 0.01 0.00 0.05 0.01 0.09 0.04 0.00 1.00 0.01
## Geology 0.00 0.00 0.01 0.84 0.01 0.95 0.10 0.08 0.00 0.00
## Total 0.00 0.00 0.00 0.50 0.00 0.69 0.22 0.00 0.00 0.00
## Regional 0.00 0.00 0.14 0.39 0.07 0.45 0.45 0.20 0.00 0.00
## Endemics 0.00 0.00 0.79 0.26 0.36 0.25 0.62 0.17 0.00 0.00
## SIE 0.00 0.00 0.02 0.39 0.01 0.50 0.01 0.01 0.00 0.00
## Regional Endemics SIE
## Area 0 0.04 0.00
## Elevation 0 0.02 0.00
## MAT 1 1.00 0.62
## MAP 1 1.00 1.00
## Temp 1 1.00 0.53
## Rain 1 1.00 1.00
## Dis 1 1.00 0.27
## Dm 1 1.00 0.53
## Geology 0 0.02 0.00
## Total 0 0.01 0.00
## Regional 0 0.00 0.00
## Endemics 0 0.00 0.01
## SIE 0 0.00 0.00
##
## To see confidence intervals of the correlations, print with the short=FALSE option
Διακρίνετε κάποιες μεταβλητές με τόσο μεγάλη συσχέτιση;
3.2.1 Αυτοματοποιημένος τρόπος
usdm::vifcor(Aegean %>% as.data.frame() %>% dplyr::select(Area:Geology), th = 0.8)## 2 variables from the 9 input variables have collinearity problem:
##
## MAT MAP
##
## After excluding the collinear variables, the linear correlation coefficients ranges between:
## min correlation ( Rain ~ Elevation ): -0.021
## max correlation ( Geology ~ Elevation ): 0.75
##
## ---------- VIFs of the remained variables --------
## Variables VIF
## 1 Area 2.6
## 2 Elevation 3.9
## 3 Temp 1.9
## 4 Rain 1.3
## 5 Dis 1.7
## 6 Dm 1.2
## 7 Geology 2.4
3.3 Οπτικοποίηση αποτελεσμάτων
Μπορούμε να οπτικοποιήσουμε τα αποτελέσματα μας και να διαπιστώσουμε εάν κάποια μεταβλητή δεν εμφανίζει κανονική κατανομή
scatterplotMatrix(Aegean[, 2:14], spread = F, smoother.args = list(lty = 2),
main = "Scatter Plot Matrix")3.4 Κανονικοποίηση των δεδομένων μας
Καθώς τα δεδομένα μας δεν έχουν κανονική κατανομή, θα πρέπει να τα κανονικοποιήσουμε. Αυτό σημαίνει ότι θα πρέπει να τα λογαριθμήσουμε. Συνεπώς, πρέπει να ελέγξουμε εάν κάποιες μεταβλητές περιέχουν το 0. Η ακόλουθη εντολή μας βοηθάει να το διαπιστώσουμε αυτό:
range(Aegean$SIE)## [1] 0 172
Κάντε το και για τις υπόλοιπες μεταβλητές (Δηλαδή:
range(Aegean$insert_IV_name_here))
3.5 Λογαρίθμηση
Εφ’όσον έχουμε ελέγξει όλες τις μεταβλητές μας, μπορούμε να τις λογαριθμήσουμε:
Aegean <- Aegean %>% mutate(Area = log10(Area), Total = log10(Total), Endemics = log10(Endemics),
Regional = log10(Regional), SIE = log10(SIE + 1))Και εάν θέλουμε να λογαριθμήσουμε όλες τις μεταβλητές που μας ενδιαφέρουν, ένας γρήγορος τρόπος είναι ο ακόλουθος:
Aegean_log <- lapply(Aegean[, 2:13], as.numeric) %>% as_tibble() %>% log10() %>%
dplyr::mutate(Island = Aegean$Island, SIE = log10(Aegean$SIE + 1)) %>% dplyr::select(Island,
everything(), SIE)3.6 Δομή των δεδομένων μας (ξανά)
Ας δούμε τα δεδομένα μας ξανά:
Aegean_log4 Κυρίως ανάλυση
4.1 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο enter
Τώρα είμαστε πλέον σε θέση να δημιουργήσουμε 4 μοντέλα πολλαπλής γραμμικής παλινδρόμησης, ένα για κάθε κατηγορία φυτικών taxa [Συνολικός αριθμός ειδών, Ενδημικά είδη, Περιφερειακά ενδημικά είδη και Τοπικά νησιωτικά είδη (Total, Endemics, Regional & SIE, αντίστοιχα)] σε σχέση με όλες τις ανεξάρτητες μεταβλητές μας2. Θα χρησιμοποιήσουμε αρχικά την μέθοδο enter:
enter1 <- lm(Total ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, data = Aegean_log)
enter2 <- lm(Endemics ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology,
data = Aegean_log)
enter3 <- lm(SIE ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, data = Aegean_log)
enter4 <- lm(Regional ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology,
data = Aegean_log)ΠΡΟΣΟΧΗ
Τα ανωτέρω μοντέλα δεν δείχνουν ποιές είναι οι μεταβλητές οι οποίες επηρεάζουν τα πρότυπα ποικιλότητας. Δείχνουν απλά τι ποσοστό εξηγείται εάν συμπεριλαμβάνονται όλες οι ανεξάρτητες μεταβλητές στο μοντέλο μας
Ας δούμε τώρα την περίληψη για κάθε ένα εξ αυτών των μοντέλων:
summary(enter1)##
## Call:
## stats::lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30188 -0.03583 -0.00143 0.04867 0.16530
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.73748 1.01870 1.71 0.094 .
## Area 0.23040 0.03857 5.97 2.3e-07 ***
## Elevation 0.01921 0.07729 0.25 0.805
## Rain 0.12354 0.11986 1.03 0.308
## Temp 0.06476 0.41385 0.16 0.876
## Dm -0.03509 0.02071 -1.69 0.096 .
## Dis 0.00981 0.02605 0.38 0.708
## Geology 0.15085 0.06651 2.27 0.028 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.088 on 51 degrees of freedom
## Multiple R-squared: 0.854, Adjusted R-squared: 0.834
## F-statistic: 42.7 on 7 and 51 DF, p-value: <2e-16
summary(enter2)##
## Call:
## stats::lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.501 -0.130 0.030 0.135 0.554
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.7112 3.0680 -0.56 0.58
## Area 0.1802 0.1161 1.55 0.13
## Elevation 0.3017 0.2328 1.30 0.20
## Rain -0.3774 0.3610 -1.05 0.30
## Temp 1.2387 1.2464 0.99 0.33
## Dm 0.0821 0.0624 1.32 0.19
## Dis -0.0205 0.0785 -0.26 0.79
## Geology 0.1331 0.2003 0.66 0.51
##
## Residual standard error: 0.26 on 51 degrees of freedom
## Multiple R-squared: 0.383, Adjusted R-squared: 0.298
## F-statistic: 4.52 on 7 and 51 DF, p-value: 0.000558
summary(enter3)##
## Call:
## stats::lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5765 -0.1688 -0.0096 0.1368 0.9108
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4846 3.4535 -0.14 0.8890
## Area 0.3987 0.1307 3.05 0.0036 **
## Elevation 0.4987 0.2620 1.90 0.0627 .
## Rain -0.1677 0.4063 -0.41 0.6816
## Temp -0.3976 1.4030 -0.28 0.7780
## Dm -0.0172 0.0702 -0.25 0.8070
## Dis 0.1056 0.0883 1.20 0.2375
## Geology -0.1697 0.2255 -0.75 0.4553
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3 on 51 degrees of freedom
## Multiple R-squared: 0.602, Adjusted R-squared: 0.547
## F-statistic: 11 on 7 and 51 DF, p-value: 2.18e-08
summary(enter4)##
## Call:
## stats::lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5901 -0.0956 0.0419 0.1268 0.4315
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.27429 2.48522 -1.32 0.194
## Area 0.23149 0.09408 2.46 0.017 *
## Elevation 0.31825 0.18856 1.69 0.098 .
## Rain 0.22590 0.29241 0.77 0.443
## Temp 1.25730 1.00964 1.25 0.219
## Dm 0.02453 0.05053 0.49 0.629
## Dis 0.00645 0.06356 0.10 0.920
## Geology 0.10430 0.16226 0.64 0.523
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.21 on 51 degrees of freedom
## Multiple R-squared: 0.557, Adjusted R-squared: 0.496
## F-statistic: 9.15 on 7 and 51 DF, p-value: 2.83e-07
4.2 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο της βηματικής πολλαπλής γραμμικής παλινδρόμησης
Με την μέθοδο αυτή θα μπορέσουμε να δούμε ποιές μεταβλητές εμπερικλείονται στο βέλτιστο μοντέλο για τον συνολικό αριθμό των ειδών
stepm1 <- stepAIC(lm(Total ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology,
data = Aegean_log), scale = 0, direction = "both", trace = 1, keep = NULL,
steps = 1000, k = 2)
summary(stepm1)Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιές μεταβλητές εξηγούν και σε τι ποσοστό τα πρότυπα της ποικιλότητας για αυτές τις κατηγορίες;
ΠΡΟΣΟΧΗ: Μας ενδιαφέρει το adjusted R squared
4.3 Ποιό μοντέλο είναι καλύτερο;
Με την βοήθεια του Akaike Information Criterion (εν συντομία: AIC), μπορούμε να διαπιστώσουμε ποιό μοντέλο είναι “καλύτερο”. Όσο μικρότερη η τιμή του AIC, τόσο καλύτερο το μοντέλο.
AIC(enter1, stepm1)Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιό από τα δύο μοντέλα για κάθε κατηγορία ειδών είναι καλύτερο;
4.4 Δημιουργία μιας λειτουργίας για τον υπολογισμό της σχετικής σημαντικότητας των μεταβλητών
relweights <- function(fit, ...) {
R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar, 2:nvar]
rxy <- R[2:nvar, 1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda^2
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta^2)
rawwgt <- lambdasq %*% beta^2
import <- (rawwgt/rsquare) * 100
import <- as.data.frame(import)
row.names(import) <- names(fit$model[2:nvar])
names(import) <- "Weights"
import <- import[order(import), 1, drop = FALSE]
dotchart(import$Weights, labels = row.names(import), xlab = "% of R-Square",
pch = 19, main = "Relative Importance of Predictor Variables", sub = paste("Total R-Square=",
round(rsquare, digits = 3)), ...)
return(import)
}4.5 Υπολογισμός της σχετικής σημαντικότητας των μεταβλητών
Μετά την περισπωμένη (~) βάζουμε όποιες μεταβλητές έχουν βγει σημαντικές από την βηματική πολλαπλή γραμμική παλινδρόμηση:
rel.1 <- lm(Total ~ Area + Dm + Geology, data = Aegean_log)
relweights(rel.1, col = "blue")| Weights | |
|---|---|
| Dm | 11 |
| Geology | 35 |
| Area | 54 |
Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιά μεταβλητή είναι σημαντικότερη για την εκάστοτε κατηγορία ειδών;
4.6 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο των ολικών υποσυνόλων
Μέχρι στιγμής, έχουμε δει δύο διαφορετικές μεθόδους, όσον αφορά την πολλαπλή γραμμική παλινδρόμηση. Όμως, δεν μπορούμε να είμαστε 100% σίγουροι ότι έχουμε φθάσει στο ορθό συμπέρασμα, καθότι ένα από τα μειονεκτήματα της βηματικής πολλαπλής παλινδρόμησης, είναι ότι μπορεί να “κολλήσει” σε ένα υποβέλτιστο πλατώ. Τον σκόπελο αυτό, μπορούμε να τον αποφύγουμε χρησιμοποιώντας την μέθοδο των ολικών υποσυνόλων:
bes.tot <- regsubsets(Total ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology,
data = Aegean_log, nbest = 4)Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).
4.7 Απεικόνιση των αποτελεσμάτων του προηγούμενου βηματος
Το γράφημα αυτό θα μας δείξει ποιός είναι ο πραγματικά βέλτιστος συνδυασμός των ανεξάρτητων μεταβλητών:
plot(bes.tot, scale = "adjr2")Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).
tot <- lm(Total ~ Area + Rain + Dm + Geology, data = Aegean_log)
summary(tot)##
## Call:
## stats::lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30141 -0.03724 0.00481 0.04501 0.16857
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9077 0.3196 5.97 1.9e-07 ***
## Area 0.2345 0.0301 7.79 2.1e-10 ***
## Rain 0.1341 0.1126 1.19 0.239
## Dm -0.0364 0.0197 -1.85 0.069 .
## Geology 0.1543 0.0636 2.43 0.019 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.085 on 54 degrees of freedom
## Multiple R-squared: 0.854, Adjusted R-squared: 0.843
## F-statistic: 78.9 on 4 and 54 DF, p-value: <2e-16
relweights(tot, col = "blue")| Weights | |
|---|---|
| Rain | 0.62 |
| Dm | 10.45 |
| Geology | 35.16 |
| Area | 53.77 |
Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιές είναι οι σημαντικότερες μεταβλητές;
4.8 Έλεγχος των αποτελεσμάτων της παλινδρόμησης
Τώρα θα ελέγξουμε κατ’αρχάς εάν οι μεταβλητές του βέλτιστου μοντέλου μας παραβιάζουν μια από τις αρχές της πολλαπλής γραμμικής παλινδρόμησης:
sqrt(vif(tot)) > 2## Area Rain Dm Geology
## FALSE FALSE FALSE FALSE
Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).
4.9 Έλεγχος πιθανοτήτων
Στο σημείο αυτό θα χρειαστεί να ελέγξουμε εάν η πιθανότητα κάθε μεταβλητής η οποία εμπερικλείεται στο τελικό μας μοντέλο είναι πραγματικά στατιστικώς σημαντική:
##--------------------------------------------------------------------------
## Let's have a look on the p-values for each of the IVs that have entered
## our final model
##--------------------------------------------------------------------------
summary(tot)##
## Call:
## stats::lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30141 -0.03724 0.00481 0.04501 0.16857
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9077 0.3196 5.97 1.9e-07 ***
## Area 0.2345 0.0301 7.79 2.1e-10 ***
## Rain 0.1341 0.1126 1.19 0.239
## Dm -0.0364 0.0197 -1.85 0.069 .
## Geology 0.1543 0.0636 2.43 0.019 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.085 on 54 degrees of freedom
## Multiple R-squared: 0.854, Adjusted R-squared: 0.843
## F-statistic: 78.9 on 4 and 54 DF, p-value: <2e-16
##--------------------------------------------------------------------------
##--------------------------------------------------------------------------
## Inside the parenthesis we include the above p-values
##--------------------------------------------------------------------------
p <- c(0, 0.239, 0.069, 0.019)
##--------------------------------------------------------------------------
##--------------------------------------------------------------------------
## Level of significance
##--------------------------------------------------------------------------
alp <- 0.05
##--------------------------------------------------------------------------
##--------------------------------------------------------------------------
## If the model's p-values are equal or smaller than the p-values we obtain
## from the following command, then they are ok
##--------------------------------------------------------------------------
p.adjust(p, method = "bonferroni", n = length(p))## [1] 0.000 0.956 0.276 0.076
##--------------------------------------------------------------------------Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).
5 Εργασία για το σπίτι
Έχετε διορία επτά (7) ημερών να στείλετε την εργασία σας σε μορφή PDF3 σε αυτό το e-mail και να απαντήσετε στα ακόλουθα ερωτήματα, αφού έχετε περιγράψει το σετ δεδομενων το οποίο έχετε χρησιμοποιήσει:
1. Πόσα νησιά έχει;
2. Πού βρίσκονται;
3. Ποιο είναι το μεγαλύτερο και ποιο το μικρότερο;
4. Ποιο είναι το ψηλότερο και ποιο το χαμηλότερο;
5. Ποιο είναι πιο κοντά στην ηπειρωτική ξηρά και ποιο είναι πιο μακρυά;
6. Ποιο έχει τα περισσότερα και ποιο τα λιγότερα είδη;
Με βάση τα αποτελέσματα σας, απαντήστε στα ακόλουθα ερωτήματα:
7. Ποιές μεταβλητές δεν εμφανίζουν κανονική κατανομή4;
8. Ποιές ενέργειες θα κάνετε προκειμένου να κανονικοποιήσετε την κατανομή των μεταβλητών αυτών;
9. Ποιές μεταβλητές εμφανίζουν μεγάλο (\(>0.8\)) συντελεστή συσχέτισης;
10. Εάν υπάρχουν κάποια ζεύγη ανεξάρτητων μεταβλητών με σ.σ. \(>0.8\), θα τις χρησιμοποιήσετε όλες; 11. Για κάθε έναν από τους τύπους μοντέλων [Enter, βηματική πολλαπλή γραμμική παλινδρόμηση και ολικά υποσύνολα], καθώς και για όλες τις κατηγορίες ειδών (Total, Endemics, SIE, Regional), να αναφέρετε το adj. r2 και τις μεταβλητές που περιέχονται στην περίληψη των μοντέλων5.
12.Συμφωνούν τα αποτελέσματα μεταξύ της βηματικής πολλαπλής γραμμικής παλινδρόμησης και της μεθόδου των ολικών υποσυνόλων; Εάν όχι, που διαφέρουν;
13. Ποιο είναι το ποσοστό συμμετοχής των μεταβλητών (για κάθε κατηγορία ειδών) στα μοντέλα της βηματικής πολλαπλής γραμμικής παλινδρόμησης και της μεθόδου των ολικών υποσυνόλων;
14. Οι τιμές των πιθανοτήτων των μεταβλητών που βρίσκονται εντός του βέλτιστου μοντέλου, είναι ίσες ή μικρότερες από αυτές του Bonferroni correction;
15. Να φτιάξετε την εξίσωση για κάθε κατηγορία ειδών, όπως αυτή προκύπτει από το βέλτιστο μοντέλο.
Τα σετ δεδομένων για την εργαστηριακή άσκηση αυτή μπορείτε να τα βρείτε εδώ.
Βλέπετε ποιά πακέτα μας ενδιαφέρουν από το επόμενο βήμα↩
Ποιές ονομάζουμε ανεξάρτητες μεταβλητές; Μπορούμε να τις χρησιμοποιήσουμε όλες; Εάν όχι, γιατί;↩
Διαφορετικά δεν θα γίνει δεκτή και δεν θα βαθμολογηθείτε↩
Να συμπεριλάβετε και το αντίστοιχο γράφημα↩
Τα αποτελέσματα να φαίνονται σε πίνακα τον οποίο θα έχετε διαμορφώσει κατάλληλα↩